The Weibull - log Weibull transition of the Interoccurrence time statistics 
in the two-dimensional Burridge-Knopoff Earthquake model 



iQ Takuma AkimotoQ and Yoji AizawJ* 



Tomohiro Hasumip Takuma AkimotoLJ and Yoji Aizawa 
Department of Applied Physics, Advanced School of Science and Engineering, Waseda University, Tokyo 169-8555, Japan 

(Dated: December 15, 2008) 

In analyzing synthetic earthquake catalogs created by a two-dimensional Burridge-Knopoff model, 
we have found that a probability distribution of the interoccurrence times, the time intervals between 
successive events, can be described clearly by the superposition of the Weibull distribution and 
the log- Weibull distribution. In addition, the interoccurrence time statistics depend on frictional 
properties and stiffness of a fault and exhibit the Weibull - log Weibull transition, which states 
that the distribution function changes from the log- Weibull regime to the Weibull regime when the 
threshold of magnitude is increased. We reinforce a new insight into this model; the model can be 
recognized as a mechanical model providing a framework of the Weibull - log Weibull transition. 

PACS numbers: 05.65.+b, 91.30.Px, 05.45.Tp, 89.75.Da 



I. INTRODUCTION 

Earthquakes are phenomena exhibiting great complexity and strong intcrmittcncy. Statistical mechanical ap- 
proaches are applied to understand the complex fault systems [l|. Although the source mechanism of earthquakes 
is still open, statistical properties of earthquakes are well described by some empirical laws, such as, the Gutenberg- 
Richter (GR) law |31 and the Omori law for aftershocks [31. The statistical properties of time intervals between succes 

nrr 

sive earthquakes, hereinafter called the interoccurrence times, have been paid much attention and discussed |5|. 
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It is shown that the 



since Bak et al. proposed the scaling law by analyzing Southern California earthquakes 
interoccurrence time distribution for the earthquakes occurring on a single fault can be described by the Weibull 
distribution, and that the Weibull exponent increases with the increase of the magnitude threshold (jj . 

Very recently, in analyzing the Japan earthquake data, we have proposed a statistical feature of interoccurrence 
times, which states that the probability distribution can be definitely written by the superposition of the Weibull 
distribution and the log- Weibull distribution . Wc have reinforced this view that the interoccurrence time statistics 
show the Weibull - log Weibull transition, which means that Weibull statistics and log- Weibull statistics coexist in the 
interoccurrence time statistics, and the dominant distribution function then changes from the log- Weibull distribution 
to the Weibull distribution as the threshold of magnitude is increased. It was found that the crossover magnitude 
from the superposition domain to the Weibull domain depends on the location on which we focused. These features 
are also observed by analyzing the Southern California and Taiwan earthquake data nj]. Meanwhile, the Weibull - log 
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FIG. 1: Schematic illustration of the 2D Burridge-Knopoff (BK) model. 



However, the theoretical background and interpretation 



Weibull transition also appeared in dynamical systems 
of this transition remain to be developed. 

Many earthquake models have been proposed and simulated to infer the source mechanism of earthquakes and 
to discuss whether the model reproduces the statistical properties of earthquakes. The Burridge-Knopoff (BK) 
model I12f is one of the theoretical models of earthquakes. Recent works on this model have been studied by many 
authors [lj. [bfl I15I lit], Il7 1 . Abaimov et al. showed that interoccurrence time distributions depended on the stiffness 
of the system by simulating the one-dimensional (ID) BK model with the static-dynamic friction law [if]]. They 
found that for a small stiffness the interoccurrence time distribution exhibits an exponential distribution, while for a 
large stiffness the interoccurrence time distribution restricted to system-wide events obeys the Weibull distribution. 
The system-wide events were defined as the events where all blocks slip during an event. In the 2D BK model, the 
present author demonstrated that the cumulative distribution of the interoccurrence time can be described by the 
power law 



which reproduces the observed behavior [18|. To the best of our knowledge, a mechanical model which 
exhibits the Weibull - log Weibull transition has not been reported. 

In this study, we attempt to understand how the Weibull - log Weibull transition and the crossover magnitude 
are influenced by a change in the major physical quantities, such as frictional properties and stiffness of a fault. In 
addition, it is worth discussing whether the BK model shows the Weibull - log Weibull transition. Thus we study 
the interoccurrence time statistics produced by the 2D BK model by changing the friction and stiffness parameters 
and the threshold of magnitude. As a result, the interoccurrence time statistics exhibit the Weibull - log Weibull 
transition. 



II. THE MODEL AND THE METHOD 



We have simulated the 2D BK model describing the relative motion of faults. As shown in fig.[TJ this model consists 
of two plates, three kinds of springs corresponding to the stress acting on the surface of faults, and blocks whose mass 
is m interconnected by the linear spring of elastic constants k^ and k v c . The blocks are connected with a fixed plate 
by the spring with spring constant k p . We assume that the slipping direction is restricted to the y-dircction. The 
equation of motion at site can be expressed by 



Vi,j = KiVi+hi + Vi-1,3 ~ 2 2/i,j) + k c(Vi,j-l + Vi.J + l - 2 Vi,j) - kpVij - F({v + iji.j]) , 



(1) 
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where yi j is the displacement, v is the plate velocity, and F is the dynamical friction force as a function of the velocity 
of the block and v. In order to rewrite this equation in a dimcnsionless form, we define the dimcnsionless position U, 
the friction function </>, and time if as 

Uij = yi.j/Do = yij/(F /kp),F(yij) = F ^(yi,j/vi),t' = uj p t = yjk p /m t, 

where F Q is the maximum friction force and V\ is a characteristic velocity. Substituting these parameters for y, t, and 
F , we obtain the non-dimensional equation of motion, 



iUu i+ i,j + Ui- ld - 2u itj ) + iKUij-! + u itj+1 - 2Uij) - Uu - <\> (27 [1/ + ir itj ] ) , (2) 



where dots indicate derivatives with respect to t' . l x (= y/k§Jkp) and l y {= \J kc/k p ) are the stiffness in the x and 
y directions, respectively, v represents the dimensionless loading velocity, which stands for the ratio of the plate 
velocity to the maximum slipping velocity . p D . ^ the dynamical friction function. As the form of 0, we use a 
velocity- weakening friction force introduced in Rcf . [19| , 

(-oo,l] U = 0, 

{l + 2 7 [C//(l-a)]} 

where 7 is a decrement in the dynamical friction force, and a is the difference between the maximum friction force 
and the dynamical friction force </>(0). To forbid a back slip, which means that a block slips in the — y direction, 
4>{U) ranges from —00 to 1, arbitrarily. This system is governed by five parameters, l x ,ly,a,v, and 7. Throughout 
this work a and v are set at 0.01, because our previous work reported that the optimal parameters of this model are 
estimated to be 1% = 1, £ = 3, a = 0.01,^ = 0.01, and 7 = 3.5 in view of the reproduction of statistical properties 
of earthquakes, for instance, the GR law with 6=1, the Zipf-Mandelbrot type power law for interoccurrence time 
statistics [17[. We calculate Eqs. © and (J3]) under the free boundary condition with the 4th order Runge-Kutta 
algorithm with the system size (N x , N y ) = (100, 25). 

In this model, a time when a block slips for the first time during an event is considered as the earthquake occurrence 
time. The nth interoccurrence time is defined as r„ = t„+i — t n , where t n and ar e the occurrence times of the 
nth and n + 1th earthquake, respectively. The interoccurrence statistics are then studied by changing l x ,ly,j, and 
the threshold magnitude m c . A seismic magnitude m in this model is defined as m = log 10 (^27 j ^ u i.j^ /1-5. Suij 
stands for the total displacement at site during an event and n is the number of slipping blocks. 

One of our goals in this work is to detect the distribution density function of the interoccurrence time P(t). For 



this purpose, we selected several kinds of distribution for P(r); the Wcibull distribution P w [9|, [20(, the log-Weibull 
distribution Pi w 0,0], the power law P pow the gamma distribution P gam (in the case of 5 = 1 in the paper Q]), 
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and the log normal distribution P;„ 2CJ, \2M which are defined by, 

CKi — 1 
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where a^Pi and /i are constants and characterize the distribution. In this time, h is fixed at 0.5. r(x) is the gamma 
function, i stands for an index number; i = 1,2,3,4, and 5 correspond to the Weibull distribution, the log-Weibull 
distribution , the power law, the gamma distribution , and the log normal distribution, respectively. Note that these 
distributions have been used as a fitting function of P(t). We will comment on the log- Weibull distribution. This 
distribution is constructed by the logarithmic modification of the Weibull distribution. Generally speaking, the tail 
of the log- Weibull distribution is much longer than that of the Weibull distribution. Previously, the log-Weibull 
distribution was derived from the chain-reaction model introduced by Huillet and Raynaud, and they then applied 
the log-Weibull distribution to fit the magnitude data in France and Japan 2l| . To maintain the statistical accuracy, 
we analyze the interoccurrence times using at least 100 events. In this work, the root mean square (rms) test and the 
Kolmogorov-Smirnov test are used in order to determine the most suitable distribution function. The rms value is 
defined as 



'E?=i(*i-^) 2 



(4) 



y n' — k 

where Xj and x\ are actual data and predicted data derived from the best-fit curve, respectively, n' is the number 
of data plots and k is the number of fitting parameters. It is well-known that the preferred distribution yields the 
smallest rms value. We calculate the rms value by use of the cumulative distribution function obtained from the 
numerical data to reduce the statistical fluctuations. Then the Kolmogorov-Smirnov test is performed in order to 
provide a more accurate confidence level, where the maximum deviation statistic Dks which is defined by 



D K s = max \ Vl - y\\ 



(5) 



where and y[ stand for the actual data of the cumulative distribution and the data estimated from the fitting 
distribution, respectively. It is well recognized that the preferred distribution has the smallest value of Dks- 



III. RESULTS AND DISCUSSION 



We attempt to trace a change in the interoccurrence time statistics by changing m c . Our previous work [17| was 
mainly focused on the no-threshold case, m c = — oo. The cumulative distributions of the interoccurrence times for 
different m c are displayed in Fig. [2] for 1% = 1,1^ = 3, a = 0.01,^ = 0.01, and 7 = 3.5. For (a), (b), and (c), m c is 
set at m c = 0.3, 0.8, and 1.1, respectively. The results of the rms value and Dks for a different distribution function 
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Interoccurrence time Interoccurrence time Interoccurrence time 

FIG. 2: The cumulative distribution of interoccurrence times in the 2D BK model for different m c . (a), (b), and (c) correspond 
respectively to the log-Weibull regime m c — 0.3, superposition regime m c — 0.8, and the Weibull regime m c — 1.1. Insets 
represent the probability density function. In (b), the solid line stands for the optimal parameter fitting by eq. (|6}, where 
ai = 1.19 ±0.005, fa = 2.21 x 10 2 ±0.56,a 2 = 7.19 ± 0.03, /3 2 = 4.31 x 10 2 ± 0.89, and p = 0.40 ± 0.04, and the symbols (x) 
are the numerical data obtained from 4609 events and 147 data points. 

TABLE I: The results of the rms value, Dks, and fitting parameters for different distribution functions. 





distribution 


Cti 


ft 


rms 


In rms 


Dks 


In Dks 




P w (i = 1) 


1.08± 0.007 


8.29xl0 2 ± 3.34 


6.4x10" 


-3 


-5.05 


0.019 


-3.96 


m c = 1.1 


P w (i = 2) 


8.05± 0.10 


1.62xl0 3 ± 13.3 


1.3x10" 


-2 


-4.34 


0.040 


-3.21 


1138 events 


Ppow — 3) 


1.56 ± 0.04 


6.36x10^ 8.42 


1.3x10" 


-1 


-2.07 


0.35 


-1.05 


70 data points 


Pgam {% = 4) 


1.02± 0.002 


8.21xl0 2 ± 4.91 


9.5x10^ 


-3 


-4.65 


0.026 


-3.65 




Pn (i = 5) 


6.33± 0.01 


0.92± 0.04 


2.4x10" 


-2 


-3.71 


0.077 


-2.56 




Pw (i = 1) 


1.31± 0.01 


5.02x10^ 0.27 


7.4x10" 


-3 


-4.90 


0.036 


-3.32 


m c = 0.3 


Pi w (i = 2) 


5.92± 0.007 


9.74x10^ 0.85 


8.5x10" 


-4 


-7.08 


0.0027 


-5.91 


19545 events 


Ppow (J> — 3) 


1.70 ± 0.04 


5.57± 0.55 


9.5x10" 


-2 


-2.35 


0.035 


-1.05 


80 data points 


Pgam {l = 4) 


1.01± 0.005 


4.86xl0 1 ± 0.90 


2.3x10" 


-12 


-3.78 


0.11 


-2.23 




Pn (i = 5) 


3.59± 0.005 


0.79± 0.01 


6.5x10" 


-3 


-5.03 


0.025 


-3.69 



for m c = 1.1 and m c = 0.3 are listed in Table U As can be seen from the table, for small m c , (e.g., m c = 0.3) 
the interoccurrence time distribution obeys the log-Weibull distribution, whereas for a large m c (e.g., m c = 1.1) 
the Weibull distribution is preferred. However, the fitting accuracy of the Weibull distribution becomes worse when 
m c is decreased. At the same time, the fitting accuracy of the log-Weibull distribution becomes worse as m c is 
increased. Hence, we think that for the intermediate case (b), the distribution can be described by the superposition 
of the Weibull distribution and the log-Weibull distribution because the Weibull components and the log-Weibull 
components of P(r) do not disappear suddenly. Actually, P(t) is well fitted by the superposition of the Weibull 
distribution and the log-Weibull distribution. Taken together, the probability distribution of the interoccurrence time 
can be expressed explicitly by the following; 



P(t) = px Weibull distribution + (1 — p) X log-Weibull distribution, 



(6) 
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FIG. 3: Change in fitting parameters as a function of m c for different 7 vaiues while fixing 1% = 1 and ly = 3 (o : 7 = 2.5, □ : 
7 = 3.5, and A : 7 = 10). The Weibull components, the log-Weibull components, and the rate of Weibull distribution are 
shown in (a), (b), and (c), respectively. 

where p means the rate of the Weibull distribution in the range, < p < 1 and depends on m c . As for p = 1, P(t) 
obeys the Weibull distribution, whereas for p = 0, P{t) follows the log-Weibull distribution. P(r) is characterized 
by five parameters, ax, 02, fa, fa, and p, while for the pure Weibull and log-Weibull regimes, the number of fitting 
parameters is two; a\ and fa for P w , and a.2 and fa for Pi w . It is reported that the interoccurrence time distribution 
can be described by the superposition of the log-Weibull distribution and the Weibull distribution in the same manner 
of cq. (J6j) by analyzing the Japan, California, and Taiwan earthquake data [9I liol]. 

A. Friction parameter 7 dependence 

Next, we focus on a change in fitting parameters by varying the friction parameter 7 and m c . For this purpose, 
we allow 7 to range from 1.0 to 10 while fixing l\ = 1 and ly = 3. Figure. [3] shows the dependence of the fitting 
parameters on m c for different 7 values (7 = 2.5 (o), 3.5 (□), and 10 (A)). As can be seen from Fig.[3](a), the Weibull 
exponent a\ gradually decreases, and the characteristic time fa increases double exponentially. Note that fa increases 
more rapidly for a large value of 7. As shown in Fig. [3] (b), the log-Weibull exponent Q2 increases linearly, and the 
characteristic time fa increases exponentially with m c . The rate of the Weibull distribution p ranges from to 1 as 
m c is increased, indicating the fact that the distribution P(t) exhibits the Weibull - log Weibull transition (see Fig. [3] 
(c)). The transition magnitude from the log-Weibull regime to the superposition regime and from the superposition 
regime to the Weibull regime is denoted respectively by m* and m** depending on 7; m* = 0.4 and m** = 1.2 for 
7 = 2.5, m* = 0.5 and m** = 1.1 for 7 = 3.5, and m* = 0.5 and m** = 0.9 for 7 = 10. We conclude that the Weibull 
- log Weibull transition holds while changing 7. 

B. Stiffness parameters 1%, and ly dependence 

In the second performance of our simulation, we focus on the relation between interoccurrence times and stiffness 
parameters, l 2 x and ly. In order to achieve this, 7 is fixed at 3.5, while l 2 x and ly are varied. As demonstrated in Fig.Q] 
(a), the change of the Weibull exponent a\ can be classified into three types; first in the case of l\ = ly = 1 (o), ati 
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FIG. 4: The relation between fitting parameters and m c for different l x and l y , whereas 7 = 3.5 (o : l x = l y = 1, □ : l x = l y = 3, 
and A : 1%. = 4 and ly = 8). The components of the Weibull distribution and the log-Weibull distribution are displayed in (a) 
and (b), respectively. In (c) the ratio of Weibull distribution is shown. 

gradually decreases as the m c is increased. Second, for l\ = ly = 3 (□), although ot\ decreases with m c in the region 
m c J$ 0.1, ai increases in the region m c gl 0.1. Finally, for l\ = 4 and ly = 8 (A), ot\ increases with the increase of m c . 
Similar tendency is observed by use of the Southern California earthquake data . f3% expands double exponentially 
as m c increases. Figure |4] (b) shows that the log- Weibull components, a-2 and P2 increase linearly and exponentially, 
respectively with m c . Finally, as clearly seen from Fig. [4] (c), the Weibull - log Weibull transition is observed because 
p changes from to 1 when m c is increased in the case, both stiffness parameters small enough, I 2 ~ 3 and l y ig 3. 
The values of the transition magnitude, to* and m** depend on l 2 x and l y ; m* = 0.4 and m** = 0.8 for l\ = l y = 1, and 
m* = —0.5 and m** = 0.6 for l\ = l 2 y = 3. It should be noted that for I 2 = 4 and ly = 8, the pure log-Weibull regime 
does not appear, because p > for any m c . In this case, the interoccurren.ee time statistics contain both the Weibull 
and the log-Weibull components, and then the dominant distribution changes from the log-Weibull distribution to the 
Weibull distribution when the threshold m c is increased. The first transition point m* cannot be determined clearly, 
but the second transition point m** is estimated to be to** = —0.1 (see Fig. |H(c)). 

C. System size dependence 

The system size dependence of the interoccurrence times is studied to discuss the finite size effect. In this time, 
the number of blocks N is changed from 625 (25 x 25) to 22500 (150 x 150), while other parameters are fixed at 
ll = l,l'y = 3, and 7 = 3.5; N is taken to be N = 625 (25 x 25), N = 2500 (50 x 50), N = 10000 (100 x 100), and 
N = 22500 (150 x 150). As clearly seen from Fig. [5l p changes from to 1 as m c is increased, suggesting the fact 
that the Weibull - log Weibull transition appears in all cases. Transition magnitudes, m* and m** are then evaluated; 
m* = 0.5 and to** = 0.7 for N = 625 (o), to* = 0.7 and to** = 1.0 for N = 2500 (A), to* = 0.5 and to** = 1.1 for 
N = 10000 (□), and to* = 0.7 and m** = 1.2 for N = 22500 (x). Thus, we conclude that the interoccurrence time 
statistics, especially the Weibull - log Weibull transition, are retained for a large system size. 
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FIG. 5: The ratio of Weibull distribution p as a function of m c for different numbers of blocks N (o : 625, A : 2500, □ : 10000, 
and x : 22500). 



D. Origin of the log- Weibull distribution 



As we mentioned, the probability distribution of interoccurrence times can be described evidently by the eq. 
However, as l\ and ly are increased, the pure log- Weibull regime becomes small (e.g., l 2 x = ly = 3) and then disappears 
(e.g., l\ = 4 and ly = 8). In this study we deduce the role of the log- Weibull distribution in view of the magnitude 
distribution. The cumulative number of earthquakes N whose magnitude is greater than or equal to m (N > m) as 
a function of m for different parameters produced by the 2D BK model are presented in Fig. [6] The arrow in Fig. [6] 
stands for the pure log- Weibull regime where p = 0. As shown this figure, when the distribution obeys the power law, 
the pure log- Weibull regime can be observed, suggesting the conjecture that the origin of the log- Weibull distribution 
in the 2D BK model is related to the power law magnitude distribution. Note that the parameter region, where the 
magnitude distribution obeys the power law globally, is limited. For l\ = l,ly = 3, and 7 = 3.5, the power law 
exponent b is b = 1.10, which is the similar to that value obtained from the earthquake data, b ~ 1. 



E. The onset mechanism of the Weibull distribution 



Next, we focus on the onset mechanism of the Weibull distribution from the viewpoint of the average event size. 
Here, the size of an event is defined as the number of slipping blocks during the event. It was shown that the 
time-interval distribution of the system- wide events obeys the Weibull distribution Ql ■ This enables us to the 
conjecture that the Weibull distribution is induced from the enhancement of the average event size, S. This conjecture 
is supported in Fig. [71 where we show the relation between the ratio p and the average event size S, and the parameter 
values correspond to the cases treated in Fig. [3] (c) . 
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FIG. 6: The cumulative number of earthquakes as a function of magnitude obtained from the 2D BK model. Circles (o), 
squares (□), triangles (A), and plus signs (+) correspond to the case of 1% = 1, ly = 3, and 7 = 3.5, Zij = 1, l y = 3, and 7 = 2.5, 
lx — 3, ly — 3, and 7 = 3.5, and l x = 4, ly = 8, and 7 = 3.5. The log-Weibull region, where p = is denoted by an arrow. All 
the plots except for the case of 1% = l,ly = 3, and 7 = 3.5 are shifted vertically for clarity. 
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FIG. 7: The relation between the average size event and the ratio of the Weibull distribution. Circles (o) and squares (□) 
correspond to the case of 1% = 1, ly = 3, and 7 = 3.5 and 1%. = l,ly = 3, and 7 = 2.5, respectively. 

IV. CONCLUDING REMARKS 

We analyzed the interoccurrence time statistics produced by the 2D BK model by varying the dynamical param- 
eters, l^,ly, and 7, for different thresholds of magnitude m c . It is found that the probability distribution of the 
interoccurrence time can be described by the superposition of the Weibull distribution and the log- Weibull distribu- 
tion. The statistics depend on 1%,, ly, and 7 and exhibit the Weibull - log Weibull transition, which states that the 
distribution function changes from the log-Weibull regime to the Weibull regime when m c is gradually increased. As 



10 



l\ and ly are increased, the log-Weibull domain becomes small and then disappears. On the contrary, the interoccur- 
rence time distribution of large magnitude events always shows the Weibull distribution. Additionally, we proposed 
a new insight into the 2D BK model; the model can be recognized as a mechanical model exhibiting the Weibull - log 
Weibull transition. In this study, it is shown for the first time that the interoccurrence time distribution exhibits the 
log-Weibull distribution, reinforcing the view that the long-range correlation hides in the 2D BK model. Thus, we 
will focus on the analysis of the spatio-temporal correlation in future. In the BK model, fault dynamics arc modeled 
as the stick-slip motion so that we infer that there is a possibility that other physical systems exhibiting the stick-slip 
motion might show the Weibull - log Weibull transition. We believe that this study provides a clue to the origin and 
the interpretation of this transition. 
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